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We study the evolution in axisymmetry of accretion disks formed self-consistently through collapse 
of magnetized hypermassive neutron stars to black holes. Such stars can arise following the merger 
of binary neutron stars. They are differentially rotating, dynamically stable, and have rest masses 
exceeding the mass limit for uniform rotation. However, hypermassive neutron stars are secularly 
unstable to collapse due to MHD-driven angular momentum transport. The rotating black hole 
which forms in this process is surrounded by a hot, massive, magnetized torus and a magnetic 
field collimated along the spin axis. This system is a candidate for the central engine of a short- 
hard gamma-ray burst (GRB). Our code integrates the coupled Einstein-Maxwell-MHD equations 
and is used to follow the collapse of magnetized hypermassive neutron star models in full general 
relativity until the spacetime settles down to a quasi-stationary state. We then employ the Cowling 
approximation, in which the spacetime is frozen, to track the subsequent evolution of the disk. 
This approximation allows us to greatly extend the disk evolutions and study the resulting outflows, 
which may be relevant to the generation of a GRB. We find that outflows are suppressed when a stiff 
equation of state is assumed for low density disk material and are sensitive to the initial magnetic 
field configuration. 
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I. INTRODUCTION 



Binary neutron star coalescence has been proposed for 
many years as an explanation of short-hard GRBs P, Q ■ 
Possible associations between short GRBs and elliptical 
galaxies reported recently @ make it unlikely that short 
GRBs are related to collapse of massive stars and super- 
novae. The merger of compact-object binaries (neutron 
star-neutron star or black hole-neutron star) is now the 
favored hypothesis for explaining short GRBs. Accord- 
ing to this scenario, the merger results in the formation 
of a stellar-mass black hole surrounded by a hot accre- 
tion torus which contains ~ 1-10% of the total mass of 
the system. Energy extracted from this system, either by 
MHD processes or neutrino-radiation, powers the GRB 
fireball. The viability of this model depends in part on 
the presence of a sufficiently massive accretion disk fol- 
lowing collapse and on whether the accretion flow pro- 
duces sufficiently energetic outflows. 

Instead of directly collapsing to a black hole, some bi- 
nary neutron star (NS) mergers could form hypermassive 
neutron stars (HMNSs) as an intermediate state. Such 
stars are supported against collapse by strong differen- 
tial rotation 4, |^, which naturally arises from the 
merger 0, i, H. The latest binary NS merger simula- 
tions in full general relativity [l^, [ll|, [l^ have confirmed 
that HMNS formation is indeed a possible outcome. We 
note, however, that HMNSs could also result from core 
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collapse of massive stars, since collapse generates strong 
differential rotation [ll,[l3 (see also [15]). 

Differentially rotating stars tend to approach rigid ro- 
tation when acted upon by processes which transport an- 
gular momentum. HMNSs, however, cannot settle down 
to rigidly rotating NSs since their masses exceed the max- 
imum allowed by rigid rotation. These objects thus un- 
dergo 'delayed' collapse to a black hole (BH). Several 
processes can act to transport angular momentum and 
drive the HMNS to collapse. For example, previous cal- 
culations in full general relativity have modeled HMNS 
evolution driven by viscous angular momentum trans- 
port |1Q] and by angular momentum loss due to gravita- 
tional radiation pdj . 

HMNS collapse due to magnetic field effects has re- 
cently been studied numerically in [13, [H, [ll] using 
codes for evolving magnetized fluids in full general rel- 
ativity [IS mi (see also [HHl,!!!). The secular angu- 
lar momentum transport in this case is provided by two 
primary MHD effects. The magnetic winding effect [1^ 
refers to the twisting of "frozen-in" magnetic field lines 
by a shear flow. The resulting magnetic stresses trans- 
port angular momentum so as to drive the fluid toward 
uniform rotation. This magnetic braking occurs on the 
Alfvcn timescale [26| , which is typically much longer than 
the dynamical time. In contrast to this smooth win ding 
process, the magnetorotational instability (MRI) [27ll28j] 
leads to exponential growth of fleld line distortions on a 
timescale comparable to the rotation period. The nonlin- 
ear outcome of this instability is MHD turbulence, which 
enhances angular momentum transport. Thus, magnetic 
braking and the MRI ultimately lead to collapse of the 
HMNS and the formation of a hot, magnetized accretion 
disk surrounding the central BH. This hot disk produces 
copious vi> pairs, and the resulting annihilation energy 
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is roughly commensurate with the requirements for a 
short-hard GRB (as long as the emission is somewhat 
beamed) 

In addition to i/P pairs, MHD processes have also been 
suggested as mechanisms of powering GRBs. In this pa- 
per, we explore the possibility of MHD-induced jet forma- 
tion by simulating the self-consistent formation of disks 
through HMNS collapse in full general relativity, con- 
tinuing the disk evolution for up to ~ 2000M after the 
collapse, and examining the physical processes at work. 
The long duration of these simulations is achieved by us- 
ing the Cowling approximation (in which the spacetime 
metric is fixed) following the phase of black hole excision 
and live evolution of the spacetime metric. (We note that 
by imposing the Cowling approximation, we cannot take 
into account changes in the metric due to the rearrange- 
ment of mass in the disk and/or accretion. However, 
we find that the change in disk mass during the Cowl- 
ing phase is ^ 3% of the total mass in all of the cases 
we study. Hence, the Cowling approximation is expected 
to be fairly accurate.) We also demonstrate our code's 
ability to handle magnetized accretion flows in stationary 
spacetimes by reproducing the results of [2§| for accre- 
tion of a magnetized Fishbone-Moncrief torus [sO] onto a 
fixed Kerr BH. The good agreement of our results with 
the published results provides confidence in our code's 
ability to handle complex MHD accretion scenarios. 

While the HMNS calculations presented here employ 
the Cowling approximation for the late phases of the evo- 
lution, we emphasize that these runs are performed using 
a code capable of handling dynamical spacetimes. Using 
the Cowling approximation allows us to probe the qua- 
sistationary behavior at much later times than presented 
in our earlier study [19]. We also note that the final 
spacetime is more general than the Kerr spacetime (our 
code allows for the presence of a massive accretion disk) . 
Accretion disk dynamics in such a spacetime have not 
been explored by previous fixed-background stationary 
spacetime GRMHD simulations. 

The evolution of accretion fiows around BHs in sta- 
tionary spacetimes and the consequent jet formation has 
been studied numerically by several groups 29, [3l|, [H, 
m, H m ill . The initial data for the studies of McKin- 
ney and Gammie [1^ and of De Villiers et al. [3l| consist 
of thick tori with weak poloidal field loops surrounding 
Kerr BHs of varying spins. The magnetic field is subject 
to the MRI, and the resulting MHD turbulence drives 
accretion onto the central BH. Magnetic field lines car- 
ried into the BH open outward and take on a stationary, 
split-monopole like structure. A relativistic, Poynting 
dominated outflow develops in this funnel region, while 
a mildly relativistic, matter dominated outflow moves 
along the outer edge of the funnel. McKinncy [34] consid- 
ered the evolution of these outflows as they propagate to 
large radius and found that the terminal Lorentz factors 
for the inner, fast jet range up to ~ 10^, easily accommo- 
dating observational constraints on GRB jets. McKinney 
also showed that the opening angle of the outflow is con- 



trolled at small radii by the corona pressure, at interme- 
diate radii by the funnel wall outflow, and at large dis- 
tances by internal magnetic stresses. In contrast, Mizuno 
et al. [35j consider a thin, Keplerian disk threaded by a 
uniform, vertical magnetic fleld. Accretion onto the BH 
again leads to a faster inner jet and a slower outer jet, 
though both are mildly relativistic. 

In this paper, we find that HMNS collapse leads to 
a magnetically dominated funnel region surrounded (in 
some cases) by a mildly relativistic, unbound outflow. 
This flow has a similar morphology to those of the mag- 
netized accretion disk simulations in stationary space- 
times described above [l^, [3l[ . Though we do not find 
a Poynting-dominated jet in the funnel, the evolution 
in this region is sensitive to the numerical handling of 
matter in hi ghly magnetically dominated regimes. As 
discussed in [23], accurate evolution is such regions is a 
significant challenge. 

Below, we first give a brief description of our formula- 
tion and numerical methods. In Section lllli we present 
results for the evolution of a magnetized torus surround- 
ing a Kerr BH, following [2^. We describe our results 
for disk evolution following HMNS collapse in Section HVl 
and summarize in Section [V] 



II. FORMULATION 

A. Basic equations and numerical methods 

The formulation and numerical scheme for our gen- 
eral relativistic, magnetohydrodynamic (GRMHD) sim- 
ulations are the same as those reported in [20|, to which 
the reader may refer for details. Here we briefly summa- 
rize the method and introduce our notation. We assume 
geometrized units {G = c = 1) except where stated ex- 
plicitly. 

We adopt the Baumgarte-Shapiro-Shibata-Nakamura 
(BSSN) formalism [s^] to evolve the spacetime metric. 
In this formalism, the evolution variables are the confor- 
mal exponent (j) = In 7/ 12, the conformal 3-metric 7^ = 
e~'*'^7ij, three auxiliary functions = — 7*-' j-, the trace 
of the extrinsic curvature if, and the tracefree part of the 
conformal extrinsic curvature Aij = e~'^'^{Kij —jijK/3). 
Here, 7 = det(7ij), and 7^ is the spatial 3-metric. The 
full spacetime metric g^i, is related to the three-metric 
Ifj-f by 7^,y = g^iv^n^n^, where the future-directed, time- 
like unit vector normal to the time slice can be written 
in terms of the lapse a and shift as n'^ = a~^(l, —/?*). 

In this paper, we assume both equatorial and axisym- 
metry, and so we only evolve the region with a; > (where 
X represents the cylindrical radius w) and z > 0. We 
adopt the Cartoon method [38| to impose axisymmetry in 
the metric evolution, and use a cylindrical grid to evolve 
the MHD and Maxwell equations. For the gauge condi- 
tions, we adopt hyperbolic driver conditions as in [39] to 
evolve the lapse a and shift 
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The fundamental variables in ideal MHD are the rest- 
mass density p, specific internal energy e, pressure P, 
four- velocity u'^, and magnetic field measured by a 
normal observer moving with a 4-velocity n'^ (note that 
B^^n^ = 0). During the evolution, we also need the three- 
velocity — /u* . The ideal MHD condition is written 
as Uf^F^^'^ = 0, where F'^'^ is the electromagnetic tensor. 
The tensor F^'^ and its dual in the ideal MHD approxi- 
mation are given by 



I 



p* 



(1) 
(2) 



where e^jya/j is the Levi-Civita tensor. Here we have in- 
troduced an auxiliary magnetic 4- vector 6^ = i?^^/-\/47r, 
where B'^^^ is the magnetic field measured by an observer 
comoving with the fluid and is related to B^^ by 



(u) 



((5f ^ + u^'u^)B'' 



The energy-momentum tensor is written as 



rp _ ^Fluid I rrnEM 



(3) 



(4) 



where Tj^""^ and denote the fluid and electromag- 
netic pieces of the stress-energy tensor. They are given 

by 



(5) 



T, 



EM 



(6) 



where h = l+ £ + P/p\s the specific enthalpy, and b^ 
b^b^. Hence, the total stress-energy tensor becomes 



T^,y = {ph 



b^)u^Ui, 



P 



9t^u- b^b^ ■ (7) 



The magnetic pressure is defined as suggested by the sec- 
ond term in the above equation: Pmag = b^/2. 

In our numerical implementation of the GRMHD and 
magnetic induction equations, we evolve the following 
conserved variables: 



P* = P^tiU^, 

B' = ^B\ 



(8) 

(9) 
(10) 

(11) 



The evolution equations are integrated in conservative 
form using a high-resolution shock-capturing (HRSC) 
scheme. Specifically, we use the monotonized central 
(MC) scheme [iOl for data reconstruction and the HLL 
(Harten, Lax and van-Leer) scheme [4l| to compute the 



flux. The magnetic field must satisfy the no monopole 
constraint diB^ = 0. Thus, we adopt the flux constrained 
transport (flux-CT) scheme In this scheme, the in- 
duction equation is differenced in such a way that a sec- 
ond order, corner-centered representation of the diver- 
gence is preserved as a numerical identity. As in [l9l | , we 
apply outer boundary conditions on the primitive vari- 
ables p,P,v\ and B^. Outflow boundary conditions are 
imposed on the hydrodynamic variables (i.e., the vari- 
ables are copied along the grid directions with the con- 
dition that the velocities be positive or zero in the outer 
grid zones). The magnetic field is linearly interpolated 
onto the boundaries. Finally, the conserved variables are 
recomputed on the boundary. 

At each timestep, the primitive variables (p, P, w') 
must be computed from the evolution variables 
(p*,-?, S'i). This is done by numerically solving the al- 
gebraic equations ©-(HU]) together with an equation of 
state (EOS), P = P{p,e). We perform evolutions with 
two types of EOS. For the Fishbone-Moncrief disk [SOj in 
Section Hill as well as one of the HMNS models (star A, 
see Section HV Ap ■ we use the F-law EOS: 



P = (F- l)p£ 



(12) 



The corresponding cold EOS is a simple polytrope, 
Pcoid = Kp^, where X is a constant. For the second 
HMNS model (star C in Section llV A|l . we assume a more 
realistic hybrid EOS [l^, HH, in which the total pressure 
is written as a sum of cold and thermal parts 



P = P 



cold 



Pt 



th 



(13) 



The cold contribution to the pressure depends only on 
the density, and is defined as follows: 



Pcold — 



A'lp^i for p < pn 
K2p^'^ for p > Pn 



(14) 



We set Fi = 1.3, F2 2.75, Ki = 5.16 x lO" cgs, K2 = 
■f^iPni^^': and pnuc = 1.8 X 10" g/cm^. This EOS has 
the desirable property that the dependence of pressure on 
density stiffens above nuclear density and the resulting 
maximum Tolman-Oppenheimer-Volkov mass (2.01Mp,) 
is in line with predictions of realistic EOSs [43|. Shock 
heating will increase the pressure above its cold value at 
a given density given by Eq. p4)) . This is reflected by 
the thermal contribution to the pressure: 

Pth = (Fth - l)p£th, (15) 
and Ecoid is the specific internal 



energy 



where £th — £ — £coid, ccow unc c 
■ consistent with the cold pressure: 

ecoid(p) = - J Pcold (p) d Q 



(16) 



In this paper, we take Fth = 1-3. 

In low-density and/or highly magnetically dominated 
regions, we find that the inversion procedure for obtain- 
ing the primitive variables (p, P, w') from Eqs. ©-(HI 
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sometimes returns an unphysical solution (i.e., a solu- 
tion with negative pressure). In such grid zones, we ap- 
ply the fix suggested by Font et al. [4J], which consists 
of replacing the energy equation (fTU]) by the cold EOS, 
P = Pcoid(p) when solving the system of equations. This 
substitution guarantees a positive pressure. In rare cases, 
this revised system also fails to give a solution and we re- 
pair the zone by averaging from nearby zones. (Averag- 
ing is not applied to the magnetic field, since this would 
introduce monopoles.) For a typical run with 500^ reso- 
lution, we find that < 10 zones require this secondary fix 
on a given timestep. 

The code used here has been tested in multiple rela- 
tivistic MHD simulations, including MHD shocks, non- 
linear MHD wave propagation, magnetized Bondi accre- 
tion, and MHD waves induced by linear gravitational 
waves [131 ■ We have also compared this code with 
the GRMHD code of Shibata and Sekiguchi ^ by 
performing simulations of the evolution of magnetized 
HMNSs [171, hj, and of magnetorotational collapse of 
stellar cores [45|. We obtain good agreement between 
these two independent codes. 



B. Diagnostics 

We keep track of the rest mass Mq and angular mo- 
mentum J on our grid by computing the following volume 
integrals: 

Mo = / p.d^x , (17) 
Jv 

J = / Suffix . (18) 
Jv 

Note that this formula for J is only valid in an axisym- 
metric spacetime [i^. The total rest mass Mq is con- 
served (baryon number conservation), and angular mo- 
mentum J is conserved in axisymmetry since gravita- 
tional radiation carries away no angular momentum. Our 
finite differencing scheme is designed to conserve Mq and 
J as a numerical identity. This is possible since the con- 
tinuity equation V^(/9u'^) = and momentum equation 
Wfj^T^i = can be written as 

dtp, + dj{p,v^) - , (19) 
dtS, + dj{a^T^,) = ^a^T^-'g^.^^ . (20) 

Taking the (^-component of the second equation, we ob- 
tain the equation for S^p in axisymmetry: 

dtS^ + dj{a^T^^)=0 . (21) 

Note that both Eqs. and pT|) are written in con- 
servative form with no source terms. This allows us to 
design a finite differencing scheme to conserve Mq and J 
as a numerical identity. 



In practice, the conservation of Mq and J will not be 
exact for three reasons. Most significantly, outflows from 
the computational grid remove both rest mass and an- 
gular momentum. Secondly, as described in Section [II Al 
the inversion from conserved to primitive variables re- 
quires on rare occasion a fix in which primitive variables 
are averaged from nearby cells. This averaging can affect 
the total rest mass and angular momentum. (In contrast, 
the fix suggested by Font et al, which is much more com- 
monly employed in our code, does not affect the rest mass 
and angular momentum since it does not change and 
Sip.) The third factor preventing strict conservation is 
the imposition of floor values for the rest mass density 
(see below). In particular, applying a density floor in- 
creases the total rest mass. We also find that the floor 
tends to increase the angular momentum as well. The 
rate of increase for these quantities can be judged from 
the early part of the simulation (before any outflow from 
the computational grid). Based on these rates, we find 
that the fractional increases in Mq and J due to the floor 
is at most ~ few x 10~^ for the entire duration of the 
runs. 

Soon after an apparent horizon forms, we excise the 
BH interior to continue the evolution p7l[39j. During the 
post-excision evolution, we compute the rest mass Mout 
and angular momentum Jout of the material outside the 
BH by computing integrals (|17p and (|18p over the volume 
outside the apparent horizon. This material includes the 
disk and corona, as well as any outbound material which 
may be on the grid. The irreducible mass Mi„ of the BH 
is given by Mi„ = yCd/l&Tr, where A is the surface area 
of the apparent horizon. Since J is conserved, we can 
compute the BH's angular momentum Jfi by 

Jh — J — Jout- (22) 

This would be exact only if the total angular momen- 
tum were strictly conserved. However, as mentioned 
above, strict conservation is broken by our atmosphere 
treatment, by our treatment of zones in which the primi- 
tive variable inversion fails, and because of outflows from 
the grid. We calculate that the angular momentum loss 
through the outer boundary is at most a few percent. 
Nevertheless, we assume J to be perfectly conserved af- 
ter excision and hence use the value of J just before ex- 
cision when computing Jh- The BH's mass Mh is then 
computed from the formula 

Mh = y'M2,. + (A/2Mi,,)2 , (23) 

which is exact for a Kerr spacetime, and is in accord 
with the formula derived using the isolated and dynami- 
cal horizon formalism [13]. 

Shortly after BH excision (hundreds of M), the space- 
time settles down to an approximately stationary state, 
and it is possible to define an (approximately) conserved 
energy: 

E^ - I ay^T\d^x . (24) 
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We can then define the fiuxes of rest mass, energy, and 
angular momentum across any closed two-dimensional 
surface 5 in a time slice: 



Kerr BH with M = 1 in Kerr-Schild coordinates [2£ 



Mo = (j> apv'd^T., 



where 



and 



E 
j 



(25) 
(26) 
(27) 

(28) 



f-ijk = n^e^^ijk is the Levi-Civita tensor associated 
with the three- metric 7^. We use the above formulae 
for calculating fluxes through the apparent horizon when 
the apparent horizon has a general shape. However, the 
disk simulations in Section Hill are performed with a fixed 
Kerr background metric for which the horizon surface 
is spherical in our adopted coordinates. In such cases, 
the surface S" is a sphere with radius r, and the fiux 
expressions reduce to 



Afo(r) = 


= dAp^v'^r^ 

J r— const 


(29) 


E{r) = 


- / dAa^T\, 

J r— const 


(30) 


j{r)^ 


/ dAa^T^^, 

J r— const 


(31) 



where dA = sin 9d9d(j). 

To determine whether a fluid particle is unbound, wc 
compute Uf. In a stationary spacetime, the value of 
Ut for a particle moving on a geodesic is conserved. If 
the particle is unbound, the radial velocity w'' > and 
—ut = l/\/l — > 1 at infinity. Hence and ut pro- 
vide an approximate criterion to determine whether a 
fluid element is unbound, provided that the fluid motion 
is predominantly ballistic and pressure and electromag- 
netic forces can be neglected. (We note that, in general, 
this condition is necessary but not sufficient even in the 
absence of external forces (isj.) 



III. MAGNETIZED DISK EVOLUTIONS 

In order to verify the ability of our code to handle 
magnetized accretion flows in stationary spacetimes, we 
perform runs to compare with published results of McK- 
inney and Gammie [2g| (see also [49]) for the evolution 
of a magnetized Fishbone-Moncrief (FM) tor us [3(J |. In 
particular, we consider the case referred to in [29| as the 
fiducial run. The spacetime metric corresponds to a fixed 



ds' 



1 ] dt^ 



+T,d9'^ + sin^ f. 
' Aar"^ sir? 9 



4r 



S 

S + ( 1 

dtdcj) 



-2a ( 1 + ^ ) sin^ 9drd<j) , 



drdt 



2r 



1 + 1 Ur^ 



sm 



(32) 



where E = r +a cos 9 and the spin parameter is chosen 
as a = 0.938. (In this section, units with M = 1 will be 
assumed.) The BH is surrounded by a FM torus specified 
by = 4.281 and inner disk radius rin = 6. The torus 
has an outer radius of 42 in the equatorial plane and the 
pressure maximum is located at r = 12. The FM solution 
provides the specific enthalpy distribution, from which 
the rest-mass density and pressure are derived assuming 
a cold polytropic EOS, P = Kp^, with F = 4/3. The 
constant K is chosen so that the maximum rest-mass 
density at i = is unity. The torus is evolved assuming 
the equation of state P = (F— l)pe (again with F = 4/3), 
in order to account for entropy generation in shocks. 

In the absence of magnetic fields and viscosity, the 
torus is in equilibrium. However, following [2^ . we add a 
small poloidal magnetic field by specifying the azimuthal 
component of the magnetic vector potential: 



A^ 



C[{p~ Pcut),0] 



(33) 



where the cutoff density pcut is chosen as 0.2pmax = 0.2. 
This form of A^ results in magnetic field loops confined 
within the torus. The proportionality constant deter- 
mines the strength of the magnetic field and is chosen so 
that max(P,iiag)/ niax(P) = 0.01. Thus, the dynamical 
equilibrium of the torus is only slightly perturbed by the 
addition of the magnetic field. 

We evolve these initial data assuming axial and equa- 
torial symmetry and using cylindrical coordinates (tu, z). 
In addition, we introduce a logarithmic radial coordinate 
transformation in order to concentrate zones near the BH 
event horizon. In particular, we take 



=(''-''0) 



ro = 0.2 



(34) 



and then w — {f/r)zu and z — {f/r)z. We perform three 
runs having 200^, 320^, and 400^ zones, with uniform res- 
olution in {w, z) and with the outer boundaries held fixed 
at Winax = ^max = 4.0. Thc resulting grid is non-uniform 
in ZD and z, and the outer boundary in the (tu, z)-plane 
is not square. At its point of closest approach to the ori- 
gin in the (w, z)-plane, the outer boundary is located at 
a radius r = 44.7. In contrast to our distorted rectan- 
gular grid, McKinney and Gammie [295 employ spherical 
polar coordinates with a logarithmic radius and 0-zones 
concentrated toward the disk plane. As in many hydro- 
dynamic simulations in astrophysics, we add a tenuous 
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t = 779.0 t = 1580 
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FIG. 1: Panels in the first row show density contours (solid black lines) and velocity vectors at selected times during the 
evolution of the magnetized FM torus. The dotted (red) contour lines enclose regions with unbound outflows (having —ut > 1 
and positive radial velocity). The second row shows the poloidal magnetic fleld structure at the corresponding times. The 
dotted blue contours correspond to (3 = P/Pmag = 0.01 and thus enclose magnetically dominated regions (such as the funnel 
region near the polar axis in the last two snapshots). In this and all subsequent density contour plots, levels are defined by 
p/pmax = io-0 3*-o o9 (j = 0-12). Likewise, for all plots of magnetic field lines, the field lines are drawn as contours of A,p 
according to = j4^,min + (^.^.max — ^.p,min)2/20 (i = 1-19), where j4^,max and A^p^min are the maximum and minimum values 
of A^. 



TABLE I: Resolution Study 



Resolution 


-Afo" E/Mo 


J/Mo 




200^ 


0.29 0.87 


1.66 


6 


320^ 


0.36 0.82 


1.45 


10 


400^ 


0.37 0.82 


1.35 


12 



The averaging period for the values of 
-Mo, E/Mo, and J /Mo in this table is 
500 < t < 1000. (Here, we assume units 
with BH mass Af — 1.) 

The ratio of the typical MRI 
wavelength at t = to the grid resolution 
in the neighborhood of the (gas) pressure 
maximum. 
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t/M 
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FIG. 2: Upper panel shows the rest-mass accretion rate (solid 
line), along with its time-averaged value 0.42 (dotted line). 
(Here, time averages are performed over the period 500 < 
t < 2000.) The middle (lower) panel shows the ratio of the 
accretion rate of total energy (angular momentum) to the rest 
mass accretion rate. The time averaged values (dotted lines) 
are E/Mo = 0.85 and J /Mo = 1.41. (Here, we assume units 
with BH mass A/ — 1.) 



"atmosphere" that covers the computational grid outside 
the star. Following [2^, this is done by imposing floors 
on both the rest-mass density and pressure as follows: 
Patm = 10-4r-3/2 and Patm = 3.3 X lO"'r~^/\ 

We first discuss the general features of the evolution. 
Figure [1] shows snapshots of isodensity contours and ve- 
locity vectors (in the top row), and of the poloidal mag- 
netic field lines (in the bottom row, drawn as contours 
of constant A^p). The snapshots at i = show the initial 
density distribution of the torus and the loops of mag- 
netic field confined within it. The MRI and magnetic 
braking cause angular momentum transport in the disk, 
leading to accretion onto the central BH and ejection of 
some material. Magnetic field lines carried into this cen- 
tral region are stretched upward by outgoing material, 
leading to the collimation of field lines in the region of 
the spin axis. This collimated structure (first seen in the 
snapshots at t = 779 in Fig. [1]) persists through the rest 
of the simulation. The panels in Fig. [T] showing the mag- 
netic field lines also have dotted contours corresponding 



to P = P/Prr 



10- 



These contours enclose the 



collimated magnetic field lines, showing that this funnel 
region is strongly magnetically dominated. 

The MRTdriven turbulence causes the magnetic field 
remaining in the disk to become highly tangled. The vi- 
olent motions of the disk lead to the ejection of material, 
especially near the outer edge of the collimated magnetic 
field region (i.e., the funnel wall, or funnel-corona inter- 
face as discussed in [1^). In the top row of panels, the 



dotted contours surround regions which have —ut > 1 
and positive radial velocity, roughly corresponding to re- 
gions of outflow. Such regions can be seen just outside 
the magnetically dominated funnel in the last two sets of 
snapshots in Fig. [T] 

Figure [2] shows the mass accretion rate Mo as a func- 
tion of time for the 400^ case, along with the ratios 
E/Mq, and J/Mq (where E and J are the total en- 
ergy and angular momentum fluxes through the horizon) . 
These quantities are calculated as integrals over the hori- 
zon surface using Eqs. (p9)) - (|3T|) . (We exclude from these 
integrals those zones which may have been affected by 
a failed primitive variable inversion. These failures tend 
to occur in the small region near the polar axis which is 
highly magnetically dominated. In particular, by exam- 
ining these failures at several times, we find that > 99% of 
the failures occur in regions where /? = P/Pmag < 10"'', 
while > 70% of the failures occur for (3 < 10"^.) The 
time averaged values (for the period 500 < t < 2000) 
are -Mq = 0.42, E/Mq = 0.85, and J/Mq = 1.41 and 
are shown as dotted lines in Fig. O These are in good 
agreement with the values given by McKinney and Gam- 
mie for the same averaging period: — A/q — 0.35, 
E/Mq = 0.87, and J/Mq = 1.46. Finally, we show in 
Table |T] that these characteristics of the accretion flow 
are consistent for all three resolutions considered. In ad- 
dition, the table gives the ratio of the typical unstable 
wavelength of the MRI to the spatial resolution. As an 
estimate of the typical MRI wavelength at i = 0, we take 
a typical value of 



i5n 



(35) 



in the equatorial plane, where — B/^/47rp is the 
(Newtonian) Alfven velocity and fl = v'''^ is the angu- 
lar velocity psi . [isj . (This expression is exact for Newto- 
nian flows with a Keplerian angular velocity profile.) Our 
highest resolution (400^) thus gives ~ 12 points across 
the MRI wavelength. We expect that the MRI is fairly 
well resolved for the two higher resolution runs (320^ and 
400^) and perhaps marginally resolved in the 200^ run. 
The rough quantitative agreement between our results 
and those of 29] is quite reasonable given the very differ- 
ent grid structures employed by the two codes and gives 
us confidence in the ability of our code to handle complex 
accretion flows. 



IV. RESULTS FOR HMNS EVOLUTION 

We now turn to magnetized HMNS collapse and the 
remnant accretion disk evolution. We describe the ini- 
tial data models and the three evolution cases in Sec- 
tion |TV3 We then describe the results and implications 
for jet formation from HMNS collapse. 
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A. Initial Data 

The evolution cases described below employ the two 
HMNS models referred to in [lH as stars A and C. (For 
the sake of consistency, those labels will be retained in 
the present paper.) 

Since star A is constructed using a F = 2 polytropic 
EOS, P = Kp^, this model may be scaled to any desired 
physical mass by adjusting the value of K [HO]- In con- 
trast, star C is constructed with the more realistic cold 
hybrid EOS given in Eq. This hybrid EOS intro- 

duces a physical density scale, and star C thus does not 
have the scale freedom enjoyed by star A. For definite- 
ness and ease of comparison between the evolution cases, 
we fix the value of the polytropic constant for star A in 
the discussion below to K = 1.37 x 10^g~^cm^s~^. This 
value was chosen so that stars A and C have the same 
rest mass (Mq = 2.95Mq). With this choice. Table HI] 
lists some properties of these stars. We note that the 
rest masses of stars A and C exceed the supramassive 
limits (i.e., the mass limits for rigid rotation) by 46% 
and 14%, respectively. These stars rotate very rapidly 
and are highly flattened due to centrifugal force. We also 
note that the parameters of star C are chosen in order 
to mimic the HMNSs formed through binary NS mergers 
with reahstic EOSs in [il|. 

As in previous papers (e.g, 0, [E Hi 113, ) i choose 
the initial rotation law u'^Uip = A^(r2c — where u'^ is 
the four-velocity, fie is the angular velocity along the 
rotational axis, and O = u'^/u* is the angular velocity. 
In the Newtonian limit, this differential rotation law be- 
comes 



The constant A has units of length and determines the 
steepness of the differential rotation. In this paper, A 
is set equal to the coordinate equatorial radius i?oq for 
star A and to 0.8i?cq for star C. The corresponding values 
of flcq/^c are shown in Table HIl (where ficq is the angular 
velocity at the equatorial surface). 

We must also specify initial conditions for the magnetic 
field. We choose to add a weak poloidal magnetic field to 
the equilibrium model by introducing a vector potential 
taking one of the following forms: 

A^ = AbW^mnxiP - Pent , 0) , (37) 

or 

A^ = Abw\jiax{p^/^ - pli^,, 0) , (38) 

with cutoffs pcut = 0.04pmax and Pcut = 0.04P,„ax- [Note 
that our previous study [l^ uses the form in Eq. (P7|).] 
As with Eq. l[33)) . these prescriptions result in poloidal 
loops of magnetic field confined within the stars. Since 
the physically realistic magnetic field configuration is un- 
known, we adopt this simple prescription as a first step. 



This is numerically convenient since confining the ini- 
tial magnetic field to the high density regions allows 
us to avoid strongly magnetically dominated regions in 
the initial data. The proportionality constant Ai, de- 
termines the initial strength of the magnetic field. We 
characterize the strength of the initial magnetic field by 
C = max(&V^')- We choose Ab such that C ^ 10"^- 
10~^. We have verified that such small initial magnetic 
fields introduce negligible violations of the Hamiltonian 
and momentum constraints in the initial data. For ex- 
ample, without the magnetic field, the normalized Hamil- 
tonian constraint violation due to discretization error is 
0.57% for Star C. For the strongest magnetic field case, 
C2 (see below), adding the magnetic field increases the 
constraint violation by a fraction 1.8 x 10~^. 

We discuss results of three evolutions, one with star A, 
and two with star C (labeled CI and C2). Case A is a con- 
tinuation of the star A run in [Toj starting from the point 
of excision {t = 2570M = 66.9Pc)- In the present paper, 
we evolve the system through the post-excision phase and 
then, for an extended period, in the Cowling approxima- 
tion. Thus, we consider the longer-timescale behavior of 
the same model presented in For case A, the mag- 
netic field at i = takes the form given in Eq. ([57]) , with 
C = 2.5 X 10-3, giving P^,^^ = 9.63 x lO^^ G. This 
run uses a constant atmosphere floor density /Oatm — 
10-^Pmax(t = 0) = 4.5 X 10"^ g/cm^. We also impose 
a pressure floor given by Patm = 5 x lO^^^Pmax- This 
pressure floor is quite small because it is determined by 
taking half of the cold polytropic pressure at the atmo- 
sphere density, and the EOS has F = 2. In practice, 
this pressure floor is rarely invoked. These values for the 
pressure and density floors are chosen as in pji] since this 
run is a continuation of the earlier run. 

The two runs involving star C differ only in their mag- 
netic field configurations. Case CI has an initial magnetic 
field of the form given by Eq. and C = 8.3 x 10-^ 

{B^^^^ = 1.81 X IQie G). Case C2, on the other hand, 
uses the form in Eq. ([55]) , which has the effect of shifting 
the initial distribution of magnetic field energy toward 
the outer layers of the star. Figure [3] shows the initial 
density contours and the magnetic field lines obtained 
from the two prescriptions for A^p . We choose the coeffi- 
cient Ab for case C2 such that C = 9.1 x 10"^ {B^^^ = 
1.84 X 10^^ G). Though this value of C is larger than our 
previous choices, we note that the average value of b^/P 
is 5 X 10"'^. Thus, the overall perturbation to the star 
is still small, and we do not observe any artifacts from 
it during the evolution. Both cases CI and C2 use an 
atmosphere prescription inspired by the disk evolutions 
in [29], in which floors on pressure and density fall off 
with radius according to: patm = 10~''pmax(i = 0)r~'^/^ 
and Patm = 1.37 X 10"*Pmax(i = 0)r~5/2^ rpj^ig prescrip- 
tion for the pressure floor gives an atmosphere which is 
significantly hotter than the cold pressure corresponding 
to Patm determined by the hybrid EOS. We find that us- 
ing a hot atmosphere with this EOS leads to fewer failures 
of the primitive variable inversion. The radial dependen- 
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FIG. 3: The left panel shows density contour lines for star C at f = 0. The middle and right panels show initial magnetic field 
lines for cases CI and C2, respectively. Contour lines have the same meaning as in Fig. [1] 



TABLE 11: Initial Models 



Model M{Mq) 


Mo (Mo) 


i?oq/M=' j/M^ r,ot/|VK|'' n^jQc" 




A 2.71 


2.95 


4.48 1.0 0.249 0.33 


38.4 


C 2.64 


2.95 


2.75 0.82 0.241 0.185 


15.5 



" The equatorial coordinate radius _Rcq normalized by the Arnowitt-Deser-Misner 
(ADM) mass. 

^ The ratio of the rotational kinetic energy to the gravitational binding energy. 
The ratio of the angular velocity at the equator to the central angular velocity. 
The initial central rotation period Pe- 



des are ad hoc, but the purpose of the atmosphere to 
stabilize the evolution without significantly affecting the 
physical outcome. Allowing the density floor to decrease 
with radius may further reduce the impact of the atmo- 
sphere on the physical behavior [5l|. Whereas case A 
begins with the post-excision phase after the HMNS has 
collapsed, cases CI and C2 are evolved from t — 0. 

All three cases are evolved with 500^ uniform spatial 
resolution in {m, z). In case A, the outer boundaries are 
located at 4.5i?oq = 20. IM = 80.3 km. For the runs 
with star C, the outer boundaries are placed at 5i?eq = 
13. 7M = 53.6 km. 



B. Results for Star A 

The presence of the small seed magnetic field in star A 
renders the star secularly unstable. Magnetic winding 
and the MRI transport angular momentum outward, 
leading to contraction of the core and expansion of the 
outer layers. The core eventually becomes radially un- 
stable to collapse, which occurs at i ~ 2535M — 66Pc- 
This leads to the formation of a BH surrounded by a 
magnetized accretion disk. (For a detailed discussion of 
the pre-collapse and early disk evolution of this model. 



see [11].) 

Our code quickly loses accuracy after the formation of 
the BH due to grid stretching on the BH throat. In order 
to prevent this, we excise the interior region surrounding 
the singularity from the grid at t = tox = 2570M. Exci- 
sion may not be necessary. Alternative methods have 
been suggested to handle the black hole spacetimc in 
the presence of matter [s^, HH . Whether this technique 
would be effective for the simulations described here de- 
serves further study. For the present, we will employ 
the standard excision technique. We note that, during 
the pre-excision evolution, the L2 norms of the Hamilto- 
nian and momentum constraints (as defined in ^lO*!) are 
satisfied to better than 1%. During the excision evolu- 
tion, the maximum values of the constraints are given by 
{n,W,My,M") = (0.5%, 0.5%, 3%, 0.75%). 

As discussed in (l9| . we are only able to perform the 
post-excision evolution accurately for ~ 400M. In order 
to consider the disk behavior on longer timescales, we 
freeze the spacetime metric after the BH has settled down 
to a quasi-stationary state. We then evolve the MHD 
equations in this fixed spacetime (i.e., the Cowling ap- 
proximation). We start the Cowling phase at t = 2997M. 
The duration of the post-excision phase (with live BSSN 
evolution) is thus 427M. (We note that the accretion 
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FIG. 4: Irreducible mass (Mirr) and the rest mass remaining 
outside the apparent horizon for case A, normalized by the 
ADM mass (M) at t = 0. Time is given in units of M (bottom 
axis) and in units of ms (top axis). By the end of the post- 
excision run, which lasts 427M, the BH has settled to a quasi- 
stationary state. 



disk thus reaches an essentially quasi-stationary state. 
To check that the lack of outflows in this case is not 
caused by our constant density floor, we also performed 
this evolution with floors which fall off with radius. In 
particular, we set Patm = 10^^Pmax(i = 0)r~^/^ and 
Pa^t„j = 10"i'*P,nax(t = 0)r~5/2. However, this did not 
make any qualitative difference in the outcome, and no 
significant outflows were observed. 

Finally, we plot accretion rates through the apparent 
horizon during the Cowling phase in Figure [D Aver- 
aging E, J, and Mq individually over the duration of 
the Cowling run, we find the ratios E/Mq = 0.81 and 
J /(MhMo) = 2.2. These are in rough agreement with 
Table 2 of [2^ , which gives results for disk evolutions with 
varying BH spin. Though their table does not give results 
for a BH with the specific spin parameter Jh/M^ = 0.8, 
(which is the estimate for the BH in case A), the nearby 
table entries suggest E/Mp ^ 0.9 and j/{MhMo) ^ 2. 
(Note that Mh = 1 in [2g|.) Roughly speaking, the 
J /{MhMo) ratio increases as the BH spin decreases be- 
cause the innermost stable circular orbit (from which 
material plunges into the BH) moves outward in radius 
as the spin decreases. This explains why accretion onto 
star A gives a higher value for J / (MiiMq) than the disk 
simulation in Section lllli for which the BH has higher 
spin Jh/Ml = 0.938. 



flow, which has little effect on the spacetime metric, is 
strongly time dependent when the Cowling approxima- 
tion is first imposed. The evolution of the disk itself is 
thus not yet quasi-stationary.) The Cowling phase of the 
evolution lasts for a further 1912M. At the beginning 
of the Cowling phase, the rotation period near the mid- 
dle of the disk is ~ 164M. Thus the post-excision and 
Cowling phases encompass ~ 14 rotation periods of the 
disk. We note that the change in the disk mass during 
the Cowling phase is about 3% of the total rest mass. 
The Cowling approximation does not take this change 
into account, but the approximation should be fairly ac- 
curate since the self-gravity of the disk is estimated to 
affect the dynamics by ~ (Afdisk/-^), and the error is 
thus ~ AMdisk/M ~ 3% (where AMdisk is the change in 
the disk mass during the Cowling phase). 

Figure [3] shows the evolution of the irreducible mass 
and the rest mass remaining outside the horizon dur- 
ing the post-excision phase. Following a period of rapid 
accretion (representing the final stages of collapse), the 
accretion slows down to a quasi-stationary rate. This 
transition occurs a,t t — tcx ^ 170M. We can estimate 
the BH mass Mh and angular momentum Jh at this time 
[see Eqs. ([22]) and and we find that Mh = 0.91 Af 

and Jh/M"^ — 0.79. Figure [5] shows snapshots dur- 
ing the post-excision and Cowling phases. We find no 
evidence for significant unbound outflows in this case 
(aside from a transient which can be seen in the first 
set of panels). Accretion from the disk takes place pri- 
marily near the equatorial plane, and very little mate- 
rial is churned up from the disk into a corona. The 



C. Results for Star C 

We now describe the evolution of the hybrid EOS mod- 
els CI and C2. The behavior of both models up to the col- 
lapse and BH formation is represented in Fig. [T) Secular 
MHD effects result in initially linear growth of |i?*'|max 
due to magnetic winding and sudden spurts of exponen- 
tial growth of |_Bf Imax due to the MRl (for a detailed 
discussion, sec [13 )• These MHD effects lead to collapse 
much later in the case of C2 than in the case of CI. This 
is due to the different magnetic field distributions in the 
two models. Since A^p is proportional to p^/"^ for C2, 
the magnetic field is stronger in the outer regions and 
weaker in the interior as compared with CI (see Fig. [3]). 
However, collapse is triggered by transporting angular 
momentum from the interior to the exterior. Since the 
interior magnetic field is weaker for case C2, this process 
takes longer and the collapse is delayed. Both models 
eventually form BHs surrounded by magnetized accretion 
disks. Constraint violations remain less than 2% during 
the pre-excision evolution for both of these models. 

For case CI, collapse occurs &t t ^ 939M = 61Pc- 
We excise the singularity ai t — t^x — 941M and freeze 
the spacetime metric at i = 1221A/. Thus, the live 
BSSN evolution with excision lasts for 280Af, and the 
maximum constraint violations during the post-excision 
phase are (H, 7W^, TW?', 7W^) = (0.5%, 4%, 8%, 3%). We 
evolve the system in the Cowling approximation for a 
further 1876Af after the post-excision evolution ends. At 
the beginning of the Cowling phase, the rotation period 
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FIG. 5: Post-excision and Cowling evolution of star A at selected times. The upper panels show density contours and velocity 
vectors and the lower panels show the poloidal magnetic field lines. The dotted (red) contour lines enclose regions with unbound 
outflows (upper panels) and regions with /p > 1 (lower panels). The heavy magenta circle centered on the origin marks the 
location of the apparent horizon. 



midway through the disk is 148Af , and thus the post- 
excision and Cowling phases together cover roughly 16 
periods of the disk. Following the collapse, the phase of 
extremely rapid accretion transitions to a slower rate at a 
time t — tcx ^ 46M. Estimating the BH mass and angu- 
lar momentum at this time [following Eqs. ([221) and 
gives Mh = 0.93M and J^/M^ = 0.71. We note that the 
collapse time found here differs from the collapse time of 
511M = 33Pc given for the corresponding model in [T9|. 
This is likely due to the sensitivity of the system as it 
becomes marginally stable. 

Snapshots of the evolution during the post-excision 
and Cowling phases are shown in Fig.[8l (Contours show- 
ing —ut — 1 are left out of this figure for the sake of 
clarity, but see Fig. [D) Several features of the evolu- 
tion are immediately apparent. As before, a low-density 
funnel region containing a coUimated magnetic field has 
formed along the axis and persists through the evolu- 
tion. In the first three snapshots (post-excision and early 
Cowling phases), there is a considerable amount of mate- 
rial and magnetic field in a low density corona above the 
disk. As material in this corona falls toward the central 
BH, the attached magnetic field lines make contact with 
the collimated magnetic field and reconnect. Material at- 



tached to the reconnected field line is then driven outward 
along the corona- funnel boundary by magnetic buoyancy. 
(This process is easiest to see in an animation, but can 
also be seen, for example, in the t — t^x = 280Af snap- 
shot of Figure [5] at x ~ 20 km and z ~ 47 km. There, an 
unbound blob of material attached to reconnected field 
lines moves toward the edge of the grid.) Though we 
have not included a physical model for resistivity in our 
code, this reconnection is ubiquitous in HRSC codes for 
ideal MHD and physical reconnection is expected to op- 
erate in systems with MHD turbulence such as this (see, 
e.g. M)- 

This process leads to an intermittent outflow along the 
corona-funnel wall boundary. The outflow is mildly rela- 
tivistic, with typical Lorentz factors ranging between 1.2 
and 1.5. In Fig. llOi we plot the approximate maximum 
asymptotic Lorentz factors associated with the outflow 
for both the excision and Cowling phases. The time- 
averaged value for the maximum is 1.2. We find that this 
outfiow dies down as the Cowling evolution proceeds (see 
Fig. [13] below). As may be seen in the last three panels of 
Fig. [51 the corona gradually empties of matter and mag- 
netic fields, leaving a quasi-stationary disk surrounding 
the BH. Without the turbulent driving of reconnection 
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FIG. 6: Fluxes through the apparent horizon during the Cowl- 
ing evolution of star A. The upper panel shows the negative of 
the rest mass flux in units of solar masses per second (where 
Mo < indicates flow into the horizon). The middle panel 
shows separately the matter and electromagnetic contribu- 
tions to the total energy flux (-Ema and -Eem respectively), 
normalized by Mo- Note that the electromagnetic contribu- 
tion has been multiplied by a factor of 100. The lower panel 
shows the electromagnetic and matter contributions to the an- 
gular momentum flux through the horizon ( Jma and Jem re- 
spectively), normalized by MhMo to yield a non-dimensional 
ratio. The electromagnetic contribution has been multiplied 
by a factor of 10. 



across the funnel-corona boundary, the outflow largely 
ceases. In addition, the draining of pressure in the corona 
allows the magnetic field in the funnel to expand outward, 
as may be clearly seen in the last panel of Fig. [S] The 
corona thus plays a role in confining the magnetically 
dominated funnel region, as found by McKinney (34| . It 
is possible that the corona would be sustained in a 3D 
evolution, for which the turbulent magnetic field growth 
is not limited by the anti-dynamo theorem [55| . With our 
current computational resources, it would not be feasible 
to perform this run in 3D with high enough resolution to 
capture the MRI. 

Case C2 undergoes collapse at t ~ 3003Af = 194Pc , 
and we excise the singularity at t — tox — 3008M. The 
Cowling phase begins at t = 3251M, giving a live BSSN 
post-excision run of duration 243M. During the post- 
excision phase, the maximum normalized Hamiltonian 
and momentum constraints are {H, , , M^) = 
(0.4%, 1, 5%, 13%, 3%). The Cowling phase continues for 
an additional 1282M beyond the end of the post-excision 
run. (A longer Cowling evolution is not necessary since 
the matter reaches an essentially quasi-stationary state.) 
The rotation period midway through the disk at the be- 
ginning of the Cowling phase is ^ 67M, and thus the 



FIG. 7: Pre-excision evolution for cases CI (solid black line) 
and C2 (red dashed line). From top to bottom, the quantities 
plotted are: (i) the central rest-mass density normalized by 
its initial value, (ii) the central lapse, (iii) the maximum value 
of \B^ \ in units of gauss (G), and (iv) the maximum value of 
IB^l in units of lO" G. Collapse occurs in case CI much 
earlier than in case C2. 



post-excision and Cowling phases together cover roughly 
22 periods of the disk. (Note that the disk in this case 
is considerably more compact than that of CI.) Follow- 
ing the collapse, the transition from extremely rapid to 
slower accretion occurs at t ~ 3055M = 197Pc- Estimat- 
ing the BH mass and angular momentum at this time 
gives Mh = 0.97M and Jh/M^ = 0.77. 

Selected snapshots from the post-excision and Cowling 
phases of evolution of case C2 are shown in Fig. [TT] In 
contrast to case CI, no significant outflows are seen for 
this case (except for a transient seen in the first panel). 
The corona region contains much less material and mag- 
netic field at the beginning of the post-excision evolution 
than in case C2. This allows the funnel region to expand, 
and it rapidly assumes the structure seen in the last two 
sets of panels. The relative absence of material in the 
corona is likely due to the much more rapid accretion 
following collapse. This can be seen in Fig. [121 which 
shows the evolution of the irreducible mass and the mass 
remaining outside the horizon in the two star C cases. For 
C2, Mout drops very rapidly after excision (solid black 
line), while the accretion takes place over a longer period 
of time for case CI. This may be understood from the 
difference in initial magnetic field configurations. For C2, 
more of the star's mass is initially threaded by magnetic 
field lines. Since the magnetic field remaining in the disk 
after collapse is thus stronger for case C2, the angular 
momentum transport is more efficient and the accretion 
is more rapid. 

The comparative lack of outflows in C2 can also be 
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FIG. 8: Post-excision and Cowling evolution phases of case CI, shown at selected times. The meanings of the lines are the 
same as in Fig [5] except that the ut = —1 contours have been left out for the sake of clarity. The first and last sets of panels 
represent the beginning of the post-excision phase and the end of the simulation, respectively. The panel at t — tcx = 280M 
marks the beginning of the Cowling phase. 
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FIG. 9: Dotted contours demarcating regions with —ut > 1, and velocity vectors. Note that the domain in these plots 
corresponds to the upper left-hand corner of the domain shown in Fig. [S] These snapshots demonstrate the funnel wall outflow, 
and the first two times correspond to the second and third sets of panels in Fig. |8] 



1 


1 1 1 1 1 


1 1 1 1 


1 1 1 


III 


1 1 1 


II 

1 1 


i) 

1 1 1 


km 

1 1 1 1 


III 

1 1 1 





500 1000 1500 2000 



(t-t^J/M 

FIG. 10: The maximum value of \ut \ on the grid as a function 
of time since excision for CI. This maximum is restricted 
to regions with outwardly directed velocity and with den- 
sity greater than 10 times the atmosphere floor level. This 
quantity gives an approximation to the maximum asymptotic 
Lorentz factor associated with the outflow. Hence, the inter- 
mittent outflow is mildly relativistic. 



seen by comparing the rest mass fluxes through an outer 
spherical shell for cases CI and C2. This is shown in 
Fig. [131 where the shell is placed at r = 12. 2M — 
47.7 km. This spherical shell cuts through the disk as well 
as the outflow region, and thus there is a negative con- 
tribution from the inflow at lower latitudes. This overall 
inflow behavior is punctuated by intermittent outflows 
in the CI case, while the flux through this shell is much 
smaller in the C2 case. These intermittent outflows even- 
tually die down as the turbulence driving the outflow 



decays. As in case C2, the lack of outflows in case A 
is likely caused by the lack of fluid and magnetic field 
in the corona region, which is in turn due to the EOS. 
Star A uses a T = 2 law for all density regimes, whereas 
the star C EOS corresponds to F = 1.3 in the low den- 
sity regions of interest in the disk and corona. The softer 
EOS for star C allows more efficient shock heating, which 
aids the ejection of material into the corona. This inter- 
pretation is corroborated by the results of Mignone and 
McKinney [56|], who found that using F = 4/3 in FM 
torus simulations leads to more vigorous turbulence and 
a thicker corona than the same case with F ~ 5/3. 

Finally we plot the various flux quantities for the two 
star C cases in Figures [U and [121 The upper panels 
show the rest mass flux through the horizon, which again 
demonstrates that, for case C2, the accretion is very rapid 
at first but then decays rapidly. The decrease in accre- 
tion rate for CI, on the other hand, is more gradual. 
Averaging over the duration of the Cowling evolution, 
we find that E/Mq = 0.90 and j/{MhMo) = 2.2 for case 
CI. For case C2, these values are instead 0.90 and 2.7. 
Cases CI and C2 produce BHs with Jh/M^ = 0.72 and 
0.78, respectively. For magnetized FM tori surrounding 
Kerr BHs with spin parameters in the range 0.7-0.8, Ta- 
ble 2 of m suggests E/Mo 0.9 and j/{MhMo) > 2.0. 
Thus, for accretion onto the central BH following col- 
lapse of star C, the values of E/AIq and J/{MhMo) are 
in rough agreement with the results of [2^ . 



V. CONCLUSIONS 

We have used a code which solves the GRMHD 
equations in dynamical spacetimes to simulate self- 
consistent disk formation and evolution following mag- 
netized HMNS collapse. Our simulations extend the re- 
sults of [l9| by following the disk evolution for a much 
longer time, as well as by considering an alternative mag- 
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FIG. 11: Post-excision and Cowling evolution phases of case C2, shown at selected times. The meanings of the lines are the 
same as in Fig [5] The first and last sets of panels mark the beginning of the post-excision phase and the end of the simulation, 
respectively. The middle panels show an intermediate time during the Cowling phase. We note that the density feature near 
the rotation axis is not outbound or unbound, and thus does not constitute a jet. 
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FIG. 12; Irreducible mass and rest mass outside the apparent 
horizon (normalized by the ADM mass at t = 0) for cases CI 
(dashed red line) and C2 (solid black line). Notice that the 
mass remaining outside the apparent horizon decreases much 
more quickly for case C2 than for CI. 



FIG. 13: Rest mass flux through a spherical shell located at 
r = 12. 2M = 47.7 km for cases CI (solid black line) and C2 
(dashed red line). This plot covers both the post-excision and 
Cowling phases. 
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netic field geometry [see Eq. ([38])]. We evolve two HMNS 
models: star A with the simple F— law EOS of the form 
P ~ (r — l)pe with r = 2, and star C with the more 
realistic hybrid EOS given by Eqs. (IT4|) and HSl). Fol- 
lowing collapse, star A quickly reaches a quasi-stationary 
accretion state, with very little matter ejected from the 
system or churned up into the corona. No significant 
outflow is observed from this system aside from a brief 
transient. This is likely due to the stiff EOS, which sup- 
presses the formation of an extended corona. Without 



significant fiuid and magnetic field above the disk, the 
funnel wall outflow mechanism does not operate. 

We considered two cases for the star C hypermas- 
sive model. Case CI develops significant outflows along 
the boundary between the corona and the magnetically 
dominated funnel region near the axis. These outflows 
are triggered by reconnection across this boundary and 
the buoyant rising of the released field lines and mat- 
ter. These outflows die down as the corona is gradually 
accreted. In contrast, case C2 develops no significant 
outflows. This model has a more extended initial mag- 
netic field than CI, and thus the newly formed BH in the 
C2 case rapidly accretes most of the material remaining 
outside the apparent horizon, leaving little material left 
to the disk and corona. Since the outfiow mechanism 
depends on interactions between the funnel and corona, 
outflows are suppressed in this case. We thus find that 
the presence of a funnel wall outfiow is sensitive both to 
the EOS and to the initial m;ignctic field configuration. 

As described in detail in [18|, the remnant disk from 
the collapse of star C may produce enough energy to 
power a short GRB through neutrino-antineutrino anni- 
hilation alone. However, we have also shown that the 
collapse results in the formation of a magnetically dom- 
inated funnel, and the subsequent disk evolution (in the 
case of CI) leads to a funnel- wall outfiow. The similar- 
ity of this morphology to previous studies of magnetized 
accretion disks [2^, [lH suggests that a Poynting domi- 
nated outflow in the funnel region may also be expected. 
That we do not see this feature may be due to the numeri- 
cal difficulty of handling magnetically dominated regions, 
especially in a GRMHD code for dynamical spacetimes 
(as discussed in (l^l), and warrants further study. In 
any case, a fully realistic evolution in this region requires 
more sophisticated treatment of microphysical processes 
(such as pair creation) and careful consideration of re- 
gions where the ideal MHD approximation may break 
down [34., If the Blandford-Znajek process [53| does 
drive a Poynting-dominated outflow in the funnel region, 
the expected luminosity is [H, Ull 

Lbz few X 10^^ a'^{B/ 10^^ Gf{M/2.QMQf erg/s . 

(39) 

This luminosity would easily satisfy the typical energy 
budget for a short GRB and would likely dominate 
the vD pair annihilation luminosity, which should be 
L^p - 10^° ergs/s [Hiil. The mildly relativistic funnel 
wall outfiow could play a role in coUimating the inner 
fast jet [U, or even in stabilizing the jet against non- 
axisymmetric instabilities ^^\. 

We also find that the magnetically dominated funnel 
region expands at late times in the simulation as the 
corona density and pressure drop. This could affect the 
opening angle at the base of a Poynting-dominated out- 
fiow. However, this emptying of the corona region may be 
due to the assumption of axisymmetry, since the corona 
depends on turbulent churning of the disk and the disk 
turbulence must decay by the anti-dynamo theorem [ssj . 
Ultimately, self-consistent 3D simulations encompassing 
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the collapse and the subsequent disk and jet evolution 
will be required. However, our 2D results should be qual- 
itatively correct before the disk turbulence begins to die 
down. 
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